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We have developed a unified finite-size scaling method for quantum phase transitions that requires no prior 
knowledge of the dynamical exponent z. During a quantum Monte Carlo simulation, the temperature is auto¬ 
matically tuned by the Robbins-Monro stochastic approximation method, being proportional to the lowest gap 
of the finite-size system. The dynamical exponent is estimated in a straightforward way from the system-size de¬ 
pendence of the temperature. As a demonstration of our novel method, the two-dimensional S = 1/2 quantum 
XY model in uniform and staggered magnetic fields is investigated in the combination of the world-line quan¬ 
tum Monte Carlo worm algorithm. In the absence of the uniform magnetic field, we obtain the fully consistent 
result with the Lorentz invariance at the quantum critical point, z = 1, i.e., the three-dimensional classical XY 
universality class. Under a finite uniform magnetic field, on the other hand, the dynamical exponent becomes 
two, and the mean-field universality with effective dimension (2 + 2) governs the quantum phase transition. 
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I. INTRODUCTION 

Recent enhancement of the computational power has en¬ 
abled us to simulate larger-scale systems with higher preci¬ 
sion than ever before. In particular, with the help of the re¬ 
cent development of simulation algorithms for strongly corre¬ 
lated quantum systems, a number of simulations have been 
performed to elucidate the novel nature of quantum phase 
transitions, in which many-body physics plays an essential 
role [1,2]. Quantum phase transitions occur at absolute zero 
temperature, triggered by quantum fluctuations. Through the 
quantum-classical mapping, a quantum phase transition in d 
dimensions, if it is of second order, can be generally de¬ 
scribed by the critical theory as the temperature-driven phase 
transition in a (d + -dimensional classical system with the 
same symmetries, where ^ is the so-called dynamical expo¬ 
nent [3, 4]. 

A world-line quantum Monte Carlo (WLQMC) method is 
one of the most powerful tools for investigating quantum crit¬ 
ical phenomena without any bias or approximation [5, 6]. A 
quantum system in d dimensions is mapped to a classical sys¬ 
tem in (d + 1) dimensions in the WLQMC method. The sys¬ 
tem length along the additional direction, the imaginary-time 
direction, is given by the inverse temperature, 13. 

When one performs a WLQMC simulation to investigate 
quantum criticality, the choice of (3 for each system size is es¬ 
sential. The reason is that the quantum critical system can be 
extremely anisotropic even if the interactions are isotropic in 
real space. While the correlation length in the real-space di¬ 
rection diverges as ^ ^ {q —9c)~^^ that in the imaginary-time 
direction does as ^ {g — 9c)~^^^ where g is the coupling 
constant that controls quantum fluctuations, Qc the quantum 
critical point, and u the critical exponent. If the dynamical 
exponent z is one, the space-time isotropy is kept aside from 
a scale factor, or the velocity of low-energy excitation. In the 
meanwhile, there are phase transitions with a dynamical ex¬ 
ponent larger than one. The Bose-Hubbard model with ran¬ 
domness that exhibits quantum criticality with ^ > 1 has been 


extensively investigated analytically [7-10], numerically [11- 
16] as well as experimentally [17, 18]. 

Let us review the renormalization group and the scaling the¬ 
ory near a quantum critical point. Consider the scale transfor¬ 
mation with a certain length scale, b. A physical quantity, 
denoted as F, is generally transformed as 

F{g - 5e, L -\r 7 = - g,), bL-\ 6^^7 

= Ly-F{LV-{g-g,),L^p-^), 

( 1 ) 

where yp is the scaling dimension of the quantity under con¬ 
sideration. In the second line of Eq. (1), we chose h = L 
and introduced F(x, y) = F(x, 1, ^). This equation has sev¬ 
eral unknown constants, gc,yF,J^,z, and the scaling function 
F{x,y) itself. 

In order to determine the constants in the finite-size scal¬ 
ing ansatz (1), one had to repeat simulations densely in the 
three-dimensional parameter space (L^P^g), and perform a 
multi-parameter finite-size scaling analysis as in Refs. 19 and 
20. It typically requires considerable computational resources 
to scan the multi-dimensional parameter space. Instead of the 
exhaustive scanning, simulations with some assumed z were 
performed in most previous studies. The consistency was 
checked after the calculation as in Ref. 21. This approach, 
however, would be awfully inefficient in the case without 
knowledge of the value of 2 : in advance. Another approach for 
z estimation was to focus on the temperature dependence of 
the correlation length, ^ 3.1 g = gc and L = oc [14]. 

After the correlation length in the thermodynamic limit was 
extrapolated at each temperature, the low-temperature asymp¬ 
totic behavior was analyzed. This two-step procedure requires 
additional computational cost, and possibly introduces some 
uncontrollable systematic error from the extrapolations, even 
if the location of quantum critical point, g^ is known. In the 
meantime, the winding numbers of the world-lines in space 
and time directions were exploited in Refs. 22 and 23. A pa¬ 
rameter, L or 13, was interpolated so that the winding num- 
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Fig 1. (Color online) Schematic pictures of virtually (a) anisotropic 
and (b) isotropic systems, where L and 13 are the system linear length 
in the real space and the imaginary time directions, respectively. The 
blue oval in each picture depicts space-time region correlated with 
the center (red cross). After the aspect-ratio optimization, the relative 
correlation lengths, /L and / P, become almost the same. 


her squared in each direction averagely took the same value. 
However, such an interpolation is again non-trivial in a multi¬ 
parameter space and multi conditions. 

One of the most effective strategies to overcome the dif¬ 
ficulty of a multi-parameter scaling is to introduce an auto¬ 
tuning technique. A number of auto-tuning techniques have 
already been used in numerical simulations in the field of sta¬ 
tistical physics. For example, the invaded cluster algorithm 
[24] or the probability-changing cluster algorithm [25] can au¬ 
tomatically locate the critical point. The Wang-Landau algo¬ 
rithm [26] enables us to directly estimate the density of states 
of a system. 

In the present paper, we employ the stochastic approxi¬ 
mation method. Recently in Ref. 27, a method to automati¬ 
cally optimize the aspect ratio of a quantum system was pro¬ 
posed for analyzing quantum criticality under strong spatial 
anisotropy. The relative correlation length, 
where and are the correlation length and the system 
size in a direction, respectively (a = x, y, or t for two- 
dimensional systems and Lr = /S), was adjusted for mak¬ 
ing the system virtually isotropic SiS Rx : Ry Rr ^ ^ ' 
1:1. Figure 1 schematically illustrates virtually anisotropic 
and isotropic one-dimensional quantum systems. In Ref. 27, 
the stochastic approximation scheme was applied to the stag¬ 
gered dimer antiferromagnetic Heisenberg model, and the uni¬ 
versality class of the quantum critical point was successfully 
identified in spite of the existence of strong finite-size correc¬ 
tions that easily lead a naive finite-size scaling analysis to an 
incorrect conclusion [28, 29]. Note that the tuning method 
using the correlation length [27] is applicable to general sys¬ 
tems, while the method based on the winding number [22, 23] 
works only for systems with U(l) symmetry. 

The aim of the present paper is to propose a unified finite- 
size scaling method based on the stochastic approximation 
technique for quantum criticality with general z. The rele¬ 
vant critical exponents including the dynamical exponent will 
be obtained simultaneously without any prior knowledge or 
assumption of the values. We will demonstrate our approach 
for the two-dimensional S = 1/2 quantum XY model in uni¬ 


form and staggered magnetic fields along z direction (in spin 
space). We will clarify that the dynamical exponent becomes 
two under a finite uniform magnetic field, while it does one in 
the absence of a uniform field. 

This paper is organized as follows: Sec. II introduces 
the scaling ansatz of the correlation lengths for space-time 
anisotropic systems, the Robbins-Monro stochastic approxi¬ 
mation method, and its convergence property. It is also dis¬ 
cussed how the stochastic approximation method is applied to 
the present finite-size scaling analysis. In Sec. Ill, the model 
considered in the present paper and the WLQMC method are 
introduced. The numerical results are shown in Sec. IV. Fi¬ 
nally our study is concluded in Sec. V. The technical details 
are reported in Appendices A and B. 

II. STOCHASTIC APPROXIMATION METHOD 
A. Conditions for realizing space-time isotropy 

As noted in Sec. I, for the system with z > 1, one should 
pay attention to the space-time aspect ratio when consider¬ 
ing the finite-size scaling analysis. In this section, we explain 
conditions to realize a virtually isotropic system during a sim¬ 
ulation. For simplicity, we assume that the model considered 
hereafter has no anisotropy in real space. Generalization to 
systems with spatial anisotropy is straightforward. 

Let us start by choosing F = ^ in Eq. (1). In this case, yp 
is one, i.e., 

L7/3). (2) 

Another choice, F = yields 

7(5 - 5c, L-\r^) = - 5c)i'/", Ly(3), (3) 

where we set yp = yr- At the quantum critical point, g = Qc, 
the correlation length in the imaginary-time direction exhibits 
the power law, oc in the limit of /3 ^ oc. By 

the definition of the dynamical exponent, one finds y^ = z. 
Then, dividing both sides of Eq. (3) by [3 yields 

7(5 - 5c, L-\p^)p = 7((5 - g,)L^I\ L^p), (4) 

where I; (a;, 5 ) = yir{x,y). 

Here, let us introduce two conditions, 

= -R, (5) 

where R is an arbitrarily chosen constant, and 

^r/P = Rr. ( 6 ) 

where R^ is another constant. Assume that conditions (5) and 
( 6 ) are both satisfied by tuning g and P in siniulatin^ systems 
with different system sizes. If this is the case, ^ and are kept 
constant even though they are different functions. Meanwhile, 
the set of arguments of ^ in Eq. (2) and that of in Eq. (4) are 
the same as each other. That is, the different functions ^ and 
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sharing the same arguments are kept constant at different 
system sizes. This means that each of the arguments should 
be constant if the functions have some reasonable monotonic¬ 
ity. The monotonicity of the scaling functions is expected to 
hold near a generic critical point and supported by our numer¬ 
ical calculation shown below. Then Eqs. (5) and ( 6 ) provide 
solutions, gc{L) and I3{L), for each system size: 

5e(i)-<7cOci:-i/" (7) 

and 

/3{L)ocL^. ( 8 ) 

Thus, the coupling constant, g, automatically converges to the 
critical point as L increases. Moreover, the dynamical expo¬ 
nent can be simultaneously estimated from the asymptotic L 
dependence of the inverse temperature, (3. 


B. Robbins-Monro stochastic approximation method 

In this section, we introduce an iteration procedure to ful¬ 
fill the conditions proposed in the previous section. Our task 
is to solve the system of nonlinear equations, ^/L = R and 
(3 = Rr, with respect to g and p for given L, R, and Rr. 
The solution cannot be obtained by standard iterative methods 
for nonlinear equations, such as the Newton-Raphson method. 
It is because ^ and have statistical errors coming from 
the Monte Carlo sampling that make the conventional meth¬ 
ods unstable. We thus employ the stochastic approximation 
method explained below. 

Let us see a concrete example of the stochastic approxima¬ 
tion. For simplicity, assume that g is already set to ^c- We 
estimate the optimal P that satisfies the relation 
The solution of this equation is denoted as /3c. First, one 
runs a short Monte Carlo simulation with a trial parame¬ 
ter and measures the correlation length, then calculates 
) = R^— • Next, one updates the parameter, /3, 

by using the Robbins-Monro type update procedure [30, 31] 

/3(”+i) =^(”) --y4(^(”)) (9) 

n 

with n = 1 and repeats the above until P^^^ converges to a 
certain value with increasing n = 2,3,4,'''. Here, p is a 
(constant) parameter that determines the gain of the feedback. 
Regardless of the choice of the gain, it is proved that P^'^^ 
converges to /3c in n ^ oc with probability one [30, 32]. 

As explained in Appendix A, the mean of the probability 
distribution of P^'^^ at the n-th step (denoted as converges 
as gn — Pc ^ where a is the derivative of A{P) at 

P = Pc, and the sign of p is chosen as the same with a. For 
ap < 1 / 2 , the variance of /3^’^^ at the n-th step (denoted as 
cr^) is evaluated as ^ 1/n^^^. For ap > 1/2, on the 
other hand, ^ s^/n, where s‘^ = a^p^ jplap — 1 ) is 
the asymptotic variance and a is the statistical error resulting 
from a Monte Carlo estimation of A{P). Here we should set 
p ^ 1 /a to minimize the variance (see the detailed discussion 
in Appendix A). By this choice of p, it is also guaranteed that 


the systematic error of P^'^^ decreases faster than the statisti¬ 
cal (standard) error. In actual simulations, one needs to per¬ 
form some (^10 at least) independent stochastic approxima¬ 
tion processes to estimate error bars of P and physical quanti¬ 
ties. The number of steps of each approximation process has 
to be large enough (> 10 ^ typically) for the systematic error 
to become negligible in comparison to the statistical error. 

The present stochastic approximation method can be ex¬ 
tended to multi-dimensional problems in a straightforward 
way. Below, we will apply the method to the quantum phase 
transition of two-dimensional S = 1/2 XY model in uniform 
and staggered magnetic fields in order to demonstrate the effi¬ 
ciency of the present approach and clarify the quantum phase 
transitions. 


III. MODEL AND QUANTUM MONTE CARLO METHOD 

A. S — 1/2 quantum XY model in uniform and staggered 
magnetic fields 


The Hamiltonian of the two-dimensional S = 1/2 quan¬ 
tum XY model in uniform and staggered magnetic fields is 
defined as follows: 




{hk) 


sh 


( 10 ) 

where (Sj ) is the -component raising (lowering) op¬ 
erator at site j, (j, k) denotes a pair of nearest-neighboring 
spins, and (h^) is the amplitude of the uniform (staggered) 
magnetic field. Here we consider the square lattice of linear 
extent L with the periodic boundary conditions, and the lat¬ 
tice is bipartite with even L. If site j belongs to one of the 
sublattices, a(j) takes zero, otherwise a{j) = 1 . 


This model can be mapped to the hard-core boson 
model with the uniform and the staggered chemical poten¬ 
tials [33]. The phase diagram of the model consists of several 
phases [34, 35]: (i) the disordered phase that corresponds to 
the insulating or pinning phase in the boson model, (ii) the 
x^-plane ferromagnetic phase with non-zero transverse mag¬ 
netization, or the compressible superfluid phase, and (iii) the 
fully-polarized phase along or the empty (fully-occupied) 
phase. We will fix to some value and change across the 
phase boundary. When is smaller than the saturation field, 
= 2,3. phase transition from the ferromagnetic phase to the 
disordered phase occurs as increases. If is larger than 
an additional phase transition from the fully-polarized 
phase to the ferromagnetic phase occurs. When = 0, 
the particle-hole symmetry holds and the phase transition is 
known to belong to the three-dimensional XY (3D-XY) uni¬ 
versality class, i.e., z = 1. Phase transitions different from 
the 3D-XY universality with ^ > 1, on the other hand, are 
expected for 7 ^ 0 [7, 34]. 
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B. World-line quantum Monte Carlo worm algorithm 

In order to simulate the system described by Hamilto¬ 
nian ( 10 ), we used the worm (directed-loop) algorithm [36- 
38] with the continuous-time path-integral representation. In 
the continuous-time representation, we introduce imaginary 
time r as 


2’ = Tre“^’^ = Tr 


exp 



( 11 ) 


where Z is the partition function. The continuous-time formu¬ 
lation was adopted because of the convenience for calculat¬ 
ing the Fourier component of the imaginary-time correlation 
function, which we will exploit to calculate Expanding 
the exponential in the r.h.s. of Eq. (11), we insert the identity, 
1 ^^) (^^1 = 1 ’ between the operators, where { 10 ^)} is 
a complete basis set of the Hilbert space. We then obtain 



Fig 2. Example of the worm scattering process. The bold line denotes 
a string of = 1/2 states (or the path of the worm head), while 
the thin line does — —1/2 states. The dotted line expresses a 
bond operator. When the worm head, Sj, (a) arrives at leg 1 of 
the bond operator, one of the following events will occur. The head 
(b) bounces back to the way which it comes from, (c) goes straight 
and gets out from leg 3, (d) jumps to leg 4, or (e) turns to leg 2. 
The transition probabilities in the scattering process are determined 
by the matrix elements (14). In the present model, event (e) never 
occurs due to the total conservation. 


Z = 


^ rP rP 

1 + ^ ^ / dn • • • / dr„ 

n 

, ( 12 ) 


E=1 


where |0n+i) = l^i)- In our WLQMC simulation, a state in 
the basis set is the direct product of the eigenstate of the local 

operator (up or down). The Hamiltonian (10) conserves to¬ 
tal of the system and thus a space-time configuration forms 
continuous lines of up spins (or down spins), i.e., the world¬ 
lines. 

One can consider the integrand in Eq. (12) as a weight 
(probability measure) of each world-line configuration. In or¬ 
der to make a simulation efficient, the second (site) term in 
Hamiltonian (10) is included in the bond term as 


i E - SI) 


(13) 


where the factor 1/4 comes from the coordination number of 
the square lattice. The matrix elements of the combined bond 
term are expressed as 


'h^/A 0 0 O' 

0 1/2 0 

0 1/2 -)i«(-l)'^«)/4 0 

0 0 0 -h'^/4 


in) 

in) 

in)- 

in) 

(14) 


A constant larger than or equal to max(|h."|/4, |/i®|/4) needs 
to be added to the diagonal elements for ensuring the non¬ 
negativity of the weights. 

In the worm algorithm, extended world-line configurations 
are introduced. The configurations to sample in the Monte 
Carlo method consist of the original world-lines and the 
world-lines with a pair of kinks, points of discontinuity. Such 
a pair is called a worm, and each of discontinuity head or tail. 
In the present spin model, the worm is represented by the pair 


of spin ladder operators, (*5^, S^) or {Sj each of which 
is defined at a space-time point. 

The whole update process of the world-line configuration 
consists of the diagonal update and the worm update [38]. 
In the former, the diagonal bond operators in the Hamilto¬ 
nian ( 10 ) are inserted into or removed from a world-line con¬ 
figuration according to the diagonal elements. In the latter, 
first a worm, i.e., a pair of the raising and lowering operators, 
is inserted at a randomly chosen space-time point, and either 
operator is chosen as the head. The order of the ladder oper¬ 
ators is uniquely determined in the case with S = 1/2 since 
the local degree of freedom is binary (up or down). The worm 
head then moves along the imaginary-time direction until it 
arrives at a bond operator. At the operator, the worm head is 
scattered and its moving direction and/or sitting site may be 
changed stochastically according to the matrix elements (14). 
In Fig. 2, an example of the worm-scattering process is illus¬ 
trated. We choose transition probabilities so as to minimize 
the bounce probability (see Fig. 2), by breaking both the de¬ 
tailed balance of each worm-scattering process and even that 
of the whole Monte Carlo dynamics [39]. This scattering pro¬ 
cess is repeated until the worm head reaches back its own tail. 
Then the head and tail destroy each other. The worm is in¬ 
serted at several times in each Monte Carlo step. 

We will investigate the phase transition between the xy- 
plane ferromagnetic phase and the disordered phase. The or¬ 
der parameter is the transverse (off-diagonal) magnetization 
m X ov y direction. Although it is non-trivial to measure off- 
diagonal correlation in a WLQMC simulation, one can effi¬ 
ciently calculate the structure factor 

= + (15) 

3,k 

the transverse susceptibility 

^=T^(E / ‘^'^1 / dT2 5+(ri)S'^:(T2)), (16) 
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and the Fourier component of the (spatial and temporal) cor¬ 
relation functions, exploiting the virtue of the worm-update 
process [36]. Here the symbols are defined as follows: (O) = 
Tr[Oe“^^]/Z and 0(r) = . The correlation 

lengths in X, y, and r directions are then calculated from the 
Fourier components by the second-moment method [40, 41]. 
The detail of the measurement of these quantities is explained 
in Appendix B. 

At the critical point, the structure factor and the susceptibil¬ 
ity exhibit the following power-law behavior: 


So{L) (X L^, 

(17) 

xiL) (X 

(18) 


respectively. Here we introduce 0 = 2 — 2/3/z/, where P is not 
the inverse temperature but the critical exponent of the order 
parameter. The exponent of the susceptibility is denoted by 7. 
Note that in Eqs. (8), (17), and (18), one can use the quantities 
evaluated at gdL), the solution of Eqs. (5) and (6) for each 
system size L, instead of the true critical point Qc, as both 
give the same exponent. 


IV. NUMERICAL RESULTS 
A. For = 0 

First we discuss the case without the uniform magnetic 
field. We performed WLQMC simulations for system sizes 
L = 8, • • • ,64, and obtained the optimal inverse tempera¬ 
ture ^{L) and staggered magnetic field hl{L) for each L by 
solving Eqs. (5) and (6) using the stochastic approximation. 
The optimal inverse temperature /3(L) ensures that the sys¬ 
tem is virtually isotropic, and the optimal staggered magnetic 
field hl{L) gives an estimate of the critical point. The tar¬ 
get relative correlation lengths, R and were chosen as 
R = Rr = 0.5, Rxy, or 0.7, where Rxy = 0.5925 is the 
approximate value of limi^^oo i{L)/L at the critical point of 
the three-dimensional classical XY model [42]. Note that a 
particular choice of R and R^ does not introduce any bias to 
the final estimates; it affects only the speed of convergence to 
the thermodynamic limit as seen below. 

The critical strength of the staggered magnetic field, , can 
be estimated based on the asymptotic form (7), i.e., 

h%{L) = hl + c{R)L-^/'', (19) 

where c{R) is a certain constant which depends only on R 
(= Rr, here). The system-size dependence of h^L) is shown 
in Fig. 3. We extrapolated the critical point from the finite-size 
data by using the ansatz Eq. (19) and assuming the same /i® for 
three different values of R. Seven parameters, /i® and (c{R), 
iy{R)) for each R, were determined by the least squares fitting. 
(Note that also iy{R) is n fitting parameter in our analysis.) In 
order to estimate the statistical error, we performed the fol¬ 
lowing bootstrap procedure. From several hundreds Robbins- 
Monro runs, the physical quantities and the parameters [(3{L) 



uR'^xy 

Fig 3. (Color online) System-size dependence of the critical stag¬ 
gered field, h%{L), for = 0. The extrapolated value is h% = 
0.99179(3), and the lines are the bootstrapped fitting curves (see the 
body). 



Fig 4. (Color online) System-size dependence of the transverse sus¬ 
ceptibility (16), the structure factor (15), and the inverse temperature 
at h%{L) for = 0. The red squares, green circles, and blue trian¬ 
gles are the data plots for = 0.5, Rxy, and 0.7, respectively. 


and hl{L)] were obtained with some statistical error bars. 
Then the data were resampled from the Gaussian distribution 
with the estimated mean and variance. The generated samples 
were fitted for large-enough system-size data that minimize 
X^/d.o.f, where the asymptotic form (19) would approximate 
the plots well (we used the data for L e [24,64] and obtained 
X^/d.o.f. « 2.0). This procedure was repeated 4000 times, 
which yielded as many fitted functions. By taking the av¬ 
erage of the function values, we obtained a whole shape of 
hl{L) that should be asymptotically accurate, which is shown 
in Fig. 3 as the solid line for each R. Finally, the bootstrap es¬ 
timation gave the critical point hi = 0.99179(3) for = 0, 
where the number in the parenthesis denotes the standard er¬ 
ror of the estimation in the last digit(s). It is consistent with 
the previous report, hi = 0.9919(4) [14], but more precise by 
an order of magnitude. 

It should be noted that in Fig. 3 the fastest convergence is 
achieved by the choice of 77 = Rxy, i.e., the leading correc- 








6 



Fig 5. (Color online) System-size dependence of the critical stag¬ 
gered field hl{L) for = 0.5. The extrapolated value is hi = 
1.21855(2). The lines were obtained by the same bootstrap approach 
with = 0 (see the body). The inset shows the detail of the extrap¬ 
olation for larger L. 

tion seems o(l/I/^/^^^) at R = Rxy- Meanwhile, the lead¬ 
ing correction is likely in the order ofl/L^/^^^ ati? = 0.5 
and 0.7, where uxy = 0.67155 [42] is the critical exponent 
of the correlation length for the 3D-XY universality class. 
These findings indicate that this phase transition belongs to 
the 3D-XY universality class. 

In Fig. 4, we present the system-size dependence of the in¬ 
verse temperature, the static structure factor, and the trans¬ 
verse susceptibility. Assuming the asymptotic forms [Eqs. ( 8 ), 
(17), and (18)], we conclude z = 0.992(6), ^jv = 1.967(6), 
and 0 = 0.968(5). These results are fully consistent with the 
scenario of the 3D-XY universality class, 7 /z^ = 1.9620(4) 
and 0 = 0.9620(4) [42], with z = 1. The final estimates 
for 2 ;, 0, and 7 /z^, which were evaluated from the asymptotic 
behavior of the different quantities, indeed satisfy the scaling 
relation 

-f/u = 0Yz ( 20 ) 

within the error bar. 


B. For = 0.5 

Next, we discuss the critical point and the critical expo¬ 
nents for the case with finite uniform magnetic field = 
0.5. In this case we used R = 0.5, 0.6, and 0.7 for L = 
8,12,16, 20, 24, 28,32,36,40,44. Following the same proce¬ 
dure with the case for = 0, we obtained hi = 1.21855(2), 
the quantum critical point. In the fitting procedure, the data 
with L G [24,44] were used (x^/d.o.f. « 0.8). 

The system-size dependence of the physical quantities is 
shown in Fig. 6. In comparison to the case with h^ = 0 shown 
in Fig. 4, larger corrections to scaling are seen, especially for 
(3{L). To cope with the strong finite-size corrections, we took 
the following procedure: Assume we have data points at sys¬ 
tem sizes L = Li, L 2 , • • •, where Li < L 2 < • • • < L^. 



Fig 6. (Color online) System-size dependence of the transverse sus¬ 
ceptibility (16), the structure factor (15), and the inverse temperature 
at hl{L) for h^ = 0.5. The red squares, green circles, and blue 
triangles are the data plots for R = 0.5, 0.6, and 0.7, respectively. 


First, we construct triads consisting of the data with three con¬ 
secutive system sizes as (Li, I/ 2 , 1/3), (I/ 2 , 7)3,1/4), • • •, and 
(I/^_ 2 ,1/n-i, 7)^), defining Lave as the average system size 
of each triad. Next, we fit each triad with a simple power 
function, i.e., y{L) = a{Lg,yQ) x with and 

6 (Lave) the fitting parameters depending on Lave- The er¬ 
ror of b{Ls,yQ) is estimated by the bootstrap method as ex¬ 
plained above. Then, we fit 6 (Lave) with a quadratic function 
of 1 /Lave 5 and extrapolate the critical exponent in the ther¬ 
modynamic limit. In our fitting procedure of the exponent, 
we assume that the fitting function should be monotonic. 

The extrapolation results are shown in Fig. 7. For the dy¬ 
namical exponent, we estimated 2 ; = 2 . 00 ( 2 ); the effective di¬ 
mension of the imaginary-time axis changes from one to two 
by the introduction of uniform magnetic field h^. As for the 
other critical exponents, 7 /z/ = 1.99(1) and 0 = 0.01(1) 
were obtained. These values coincide with the mean-field ex¬ 
ponents, i.e., 7 /z^ = 2 and ^ = 0. This is consistent with 
the result for the dynamical exponent, ^ = 2 , by which the 
effective dimension of the critical theory becomes four, the 
upper critical dimension. Thus we have demonstrated that our 
finite-size scaling method enables us to extract the dynami¬ 
cal exponent successfully without any prior knowledge of the 
value of 2 ;. The universality of the quantum phase transition 
without particle-hole symmetry belongs to the mean-field uni¬ 
versality class. This is consistent with the discussion on the 
Bose-Hubbard model in Ref. 7. 


V. SUMMARY AND DISCUSSIONS 

In the present paper, we have presented the unified finite- 
size scaling method that works well regardless of the value 
of the dynamical exponent. During the WFQMC simulation, 
the system size in the imaginary-time direction in the path- 
integral representation is adjusted automatically so as to sat¬ 
isfy the conditions, ^/L = R and = Rr^ based on the 
Robbins-Monro stochastic approximation. This auto-tuning 
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hard-core boson system. The Lorentz invariance, 2 ; = 1, re¬ 
flects the particle-hole symmetry at half-fllling. In the case 
with = 0.5, we have concluded that the dynamical expo¬ 
nent changes to two and the other exponents take the mean- 
fleld values. This result of the mean-fleld universality is con¬ 
sistently explained by the conclusion that the dimension of the 
effective held theory is four; d -f ^ = 4, the upper critical di¬ 
mension. Our conclusion z = 2 agrees on the discussion in 
Ref. 7, in which the authors claimed that the fourth-order term 
in the effective action becomes irrelevant. 

The method proposed in the present paper is applicable also 
to models with randomness. For example, the method will be 
effective for systems whose dynamical exponent depends on 
the parameters, such as the random Ising model in random 
transverse held [20, 43], where the dynamical exponent may 
take an irrational value or even becomes infinite [44]. It is thus 
extremely difficult to analyze the properties of quantum crit¬ 
icality by the conventional strategies. By using our method, 
one does not need any assumption about the dynamical ex¬ 
ponent, which should be quite effective for the systematic in¬ 
vestigation of the randomness-driven quantum critical point 
causing extreme space-time anisotropy. 
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Fig 7. (Color online) Finite-size corrections of the critical exponents: 
(a) the dynamical exponent, z, (b) the exponent of the susceptibility, 
7 /z/, and (c) that of the structure factor, 0. Each data point represents 
the exponent obtained by the fit for each triad (see the body). 


procedure guarantees that the coupling constant converges to 
the critical point and the inverse temperature is proportional 
to for large enough L. 

We then applied the method to the two-dimensional S = 
1/2 quantum XY model in uniform and staggered magnetic 
fields. The correlation lengths were measured by the worm 
algorithm based on the continuous imaginary-time represen¬ 
tation. In the absence of the uniform magnetic held, = 0, 
our numerical results are consistent with the 3D-XY univer¬ 
sality class. This system can be mapped into the half-filled 


Appendix A: Convergence by the Robbins-Monro algorithm 

In this appendix, the time evolution of the probability 
distribution function driven by the Robbins-Monro iteration 
[Eq. (9)] is discussed. Let us consider the situation in which 
the physical quantity A{/3) is obtained by Monte Carlo sim¬ 
ulation and the distribution function of A{f3) is given by a 
Gaussian (normal) distribution written as 

1 

= / exp 

V2Tra^ 

Also, we assume that /(/3) can be expanded as 

m « a(/3 - /?c), 


(A-my 

2cr2 


(Al) 


(A2) 
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Fig 8. (Color online) Convergence of the (a) bias and (b) variance 
for solving /(/3) = exp(^) — 1 = 0 with the initial condition 
= /3init = 1. At the n-th RMS, a value (corresponding to the 
relative correlation length in the quantum-phase-transition analyses) 
is generated from the normal distribution, Note 

that a = 1 in this case. The results for p = 0.1, 0.3, 0.5, 1.0, 2.0, 
and 5.0 are shown (calculated from 10® Robbins-Monro processes). 
As shown in the text, the variance is minimized by the choice of 
p = Xja — lat large enough RMSs, while the bias convergence is 
accelerated monotonically with respect to p. 


near /3c, the zero of /(/3). Then, the asymptotic recursion 
relation between the distribution functions of and 
is written as 


n-\-l 




X exp • 


1 


!?:(/3(n) _^(n+l)) _^(/3(n) 

P 


1 2 ' 


(A3) 


In our procedure, we start from an initial condition 
— /3init). Here we assume that the distribu¬ 
tion function of ) will be approximated by a Gaussian 

for large enough n. Then the recursion relations for the mean. 


jin, and the variance, cr^, of Pn{/3^^^) are obtained as 

= (l - ^) + ^/3c (A4) 

V n J n 

respectively. Eq. (A4) can be rewritten as 

Mn+l - Pc = (l - —) (Mn - Pc) ■ (A6) 

V n J 


The absolute value of (1 — ap/n) in Eq. (A6) is less than 1 for 
sufficiently large n. Then ^ ^^c for ^ ^ co. Similarly, it 
can be proved that ^ 0 forn ^ oc. 

Next, let us assume the leading term of as , where 

as is an unknown constant and is the asymptotic variance 
ofcr^. Erom the approximation (n+l)“^^ n~^^{l—as/n) 

forn 1, the lowest-order terms in Eq. (A5) are evaluated as 

_ s^as ^ _ 2s^ap 

rj^Ots + l ^2 rj^Ots + l 


When < 1, dominates over 1/n^. Then ag = 

2ap and 




J^QiOjP 

When as = 1, on the other hand, 
1 cr^p^ 


for ap < -. 


s 

n 


n 2ap — 1 


for ap > -. 


(A8) 


(A9) 


There is no solution under the assumption as for the case 
where ap = 1/2, but it is expected that only some correction 
from ^ 1/n will appear. 

Similar discussion holds for the mean of distribution, jin- 
Assuming — /3c = with some constants k and am, 

we obtain 



n 

(AlO) 


from Eq. (A4). This results in am = ccp, and thus we have 

l^n- Pc - \z for ap > 0. (All) 

nap 

According to the asymptotic forms (A8), (A9), and (All), 
let us discuss the dependence of the final error on the num¬ 
ber of Monte Carlo steps. It is obvious from Eqs. (A8) and 
(A9) that it is better to set the gain |p| large enough so that 
ap > 1/2 is satisfied. Otherwise the convergence of the 
variance becomes slower. We will call one iteration of the 
Robbins-Monro feedback process [Eq. (9)] a Robbins-Monro 
step (RMS), and suppose that a whole calculation consists of 
A^r RMSs and each RMS has A^m Monte Carlo updates. The 
total computational cost is proportional to A^tot = 3 Vr x 3Vm • 
Erom Eq. (A9), the variance of the estimate after A^r RMSs 
is given by ctatr ~ 5^/A^r ^ ct^/A^r for ap > 1/2. Using 
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cr^ ~ Smc/^m, where is the asymptotic variance of the 
Monte Carlo estimation, we obtain 




A^tot' 


(A12) 


This means that the asymptotic variance depends only on the 
total number of Monte Carlo updates A^tot • 

Next, let us discuss the optimal choice of the gain p. For 
ap > 1/2, the convergence of pn is faster than that of 
Thus, we should minimize in Eq. (A9); then we derive 


p=-. (A13) 

a 

Fig. 8 shows the p dependence of pn and calculated from 
an exemplary case: 


/(^) = exp(/3) - 1, (A14) 

where /3c = 0 and a = 1. While the convergence of the 
mean becomes faster monotonically as p increases, the vari¬ 
ance takes a minimum value atp=l/a = l and increases 
again for larger p. In the simulation presented in the main 
text, we performed short preparatory calculations for small 
systems in order to roughly estimate a ^ a*, and set p = 1/a* 
for succeeding long runs. We performed several hundreds 
of Robbins-Monro iterations, each of which has A^m = 500 
worm updates. 


Appendix B: Measurement of off-diagonal correlation in the 
worm algorithm 

We explain the way to measure off-diagonal correlation in 
the worm algorithm. Let us begin with measuring the static 
structure factor defined in Eq. (15). This quantity is easily 
rewritten by the spin raising and lowering operators. We then 
need to evaluate the thermal average. 


respectively. In other words, Eq. (Bl) can be read as the fre¬ 
quency of events that the worm head visits the same imaginary 
time with the tail. We thus can measure the static structure 
factor by simply counting the frequency of extended config¬ 
urations that contribute to the numerator in Eq .(Bl) during 
each worm update. 

We have also evaluated the dynamic structure factor at 
imaginary frequency iuj. It is given by the canonical corre¬ 
lation function as 

C{q, iuj) = / ^^2 Sf {ti)S^ (ra) 

X exp {-i [w(t2 - Ti) + q- (n - fj)]} ), (B2) 

where iuj = 27ri/P is the lowest Matsubara frequency in our 
simulation. The spatial phase factor can be cal¬ 

culated and stored in advance. Then, when the head moves in 
the worm-update process, the imaginary-time integral is eval¬ 
uated. In simulations, every time the head reaches a bond 
operator, a part of the imaginary-time integral is performed. 
Since the r-dependent part of the integrand is simply given by 
we can evaluate the part of the integral exactly at each 
head move. The spatial phase factor is multiplied (if neces¬ 
sary). The matrix elements of the head and tail also need to 
be considered (they are simply one in the case with S = 1/2). 
The contribution to the integral at each head move is summed 
up until the head returns back to its tail. The average value of 
the summed integral for each worm insertion will provide the 
target quantity (B2). 

The transverse susceptibility (16) is expressed as x = 
(7(0,0). The correlation length can be estimated by the 
second-moment method [40, 41]; the correlation length in the 
X direction, is expressed as 


1 / C(go,0) 

\5qx\ y C{qo+5qx,Qi) 


(B3) 


(Es/s^) 


Tr 




Tre-/5^ 


(Bl) 


The numerator and denominator in Eq. (Bl) correspond to the 
partition function of extended world-line configurations with 
a worm and those of the original world-line configurations. 


where = {21^ jL^ 0) and go = (0? 0)- Similarly, the corre¬ 
lation length in the imaginary-time direction, as 


^ 1 I C{go, 0) 
f^y C{qo,iio) 


(B4) 
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